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SECTION  I 


INTRODUCTION 

Statistical  scatter  factors  to  be  used  in  the  reliability 

assessment  of  aircraft  structures  were  evaluated.  Two  programs 

were  written  for  this  purpose.  The  first  program  (A7701) 

calculates  the  expected  value  of  the  reliability  E[R']  (Eq.22), 

the  ratio  E[R']/R  where  R is  the  assigned  value  of  reliability, 

the  coefficient  of  variation  VR,  = aR,/E[R']  where  the  aRl  = 

2 2 P 

{ E [ ( R ' ) ] - (E[R'])  } 2 is  given  in  Eg.  23,  and  the  corresponding 

scatter  factors  from  Eg.  6,  S based  on  R and  SI  based  on  E[R'] . 

The  eguation  numbers  correspond  to  eguations  in  Volume  I.  The 

calculations  in  the  program  A7701  were  based  on  the  criteria 

that  the  shape  parameter  in  the  Weibull  distribution  is  known 

and  the  scale  parameter  is  found  by  maximum  likelihood  estimation 

(MLE) . In  the  second  program  (A7720)  the  scatter  factors  are 

estimated  for  the  case  where  both  the  scale  and  the  shape 

parameters  in  the  Weibull  distribution  are  unknown.  For  this 

purpose  the  Monte  Carlo  simulation  technique  is  introduced. 

^ ^ * 

The  distributions  of  parameters  u , Vq  and  Z defined  in 
Eqs . 32,  33,  36  respectively,  are  empirically  estimated.  This 
program  also  calculates  the  two  dimensional  frequencies  of  u 
and  Vq  (Eqs.  28  and  29),  the  reliability  ratio  E[R"]/R,  and 
the  coefficient  of  variation  VR„  = cjrh/E[r"]  (ctRII  = E[(R")Z]- 
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(E[R"])2}2)  where  E[R"]and  E[(R")2]  are  obtained  from  Eq.  41. 
The  following  sections  of  this  document  discuss  these  two 
programs,  the  input  data,  and  the  output. 


SECTION  II 

SCATTER  FACTOR  WITH  KNOWN  SHAPE  PARAMETER  (Program  A7701) 

A simplified  flow  chart  of  Program  A7701  is  shown  in 
Figure  1. 

2 . 1 INPUT  DATA 

There  are  no  input  data  cards  used  for  this  program. 

The  necessary  data  parameters  need  to  be  specified  at  the 
beginning  of  the  program  through  DATA  statements.  The  input 
parameters  are  as  follows:  AM  = grid  size  indicating  the  fleet 
size  of  aircraft;  A = shape  parameter  in  the  Weibull  distri- 
bution; R = fleet  reliability;  NPT  = number  of  division  points 
on  the  distribution  of  Z;  DZ  = increment  size  in  the  distri- 
bution . 

2.2  PROGRAM  DESCRIPTION 

Program  A7701  consists  of  a main  program  unit  and  two 
subroutines:  GAMMA  and  SIMPN.  The  program  listing  is  given 
in  Appendix  A. 
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SUBROUTINES 


START 

— 1 

INPUT  DATA 

J 

GAMMA  FUNCTION 

I 

R 

1 

INTEGRATE  EQ.  22 
E[R'  ] 

I 

INTEGRATE  EQ.  23 
E[(R')2] 


STANDARD  DEVIATION 
of  R' 

I 

COEFFICENT  OF 
VARIATION  OF  R' 

I 

SCATTER  FACTOR 
BASED  ON  R 

J 

SCATTER  FACTOR 
BASED  ON  E [R ' ] 


, 1. 

I i 

■I GAMMA  I 


SIMPN 


^ SIMPN  | 

I I 


END 


Figure  1.  Simplified  Flow  Chart  of  Program  A7701 
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2.2.1  The  Main  Program  Unit  (A7701) 

The  main  program  unit  calculates  the  expected  value  of 

E[  R']  of  R'  from  Eq.  22,  the  ratio  E[R']  /R  where  R is  given 

as  the  input  parameter,  the  variance  VAR[R']  from  Eqs . 22  and 

2 ~ o 

23  as  VAR  [ R ’ ] = ( aRl  ) =E  [ (R'  ) 1 ] - (E  [R'  ] ) * and  the  scatter  factors  , 
S based  on  R and  SI  based  on  E[R']  . The  values  for  the 
Gamma  function  and  the  numerical  integration  necessary  to 
complete  Eqs.  22  and  23  are  obtained  from  subroutines 
GAMMA  and  SIMPN. 

2.2.2  Subroutine  GAMMA  (N,  GA) 

Subroutine  GAMMA  evaluates  the  Gamma  function 

r(n)  = (n-1) ! for  specified  positive  integer  values  of  n. 

The  arguments  for  this  subroutine  are: 

N = n,  argument  of  the  Gamma  function 
GA  = output  value  of  the  Gamma  function 

2.2.3  Subroutine  SIMPN  (AF,  NPOINT,  DSTEP , AREA) 
Subroutine  SIMPN  is  a numerical  integration  routine  which 

utilizes  Simpson's  rule  to  obtain  area  for  a specified  function. 

The  arguments  for  this  subroutine  are: 

AF  = function  to  be  integrated 

NPOINT  = number  of  points  or  subdivisions  on  the 
function  AF 

DSTEP  = integration  step 
AREA  = area  from  the  integral 
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2.3  PROGRAM  OUTPUT 


Program  A7701  prints  out  the  output  as  shown  in  Table  1. 


TABLE  1.  Results  from  Program  A7701 


sample  SIZE  N = 1 

WEIBULL  SHAPE  A<  1)=  .50000  INCREMENT  ON  Z DZ=  .02500 

FLEET  RELIABILITY  LEVEL  R=  .50000 


E(R")  = 

• 

42468E+00 

E ( K* 

')/R=  . 84  93  7E  + 0 0 

VAR (R") = 

• 

66335E-0 1 

STD 

= .25756E+00 

VR"  = 

•60647E+00 
FLEET  SIZE 

S.F. 

BASED  ON  R 

S.F. 

BASED  ON  E(R") 

M = 

1. 

5 = 

. 10000E  + 01 

SI  = 

• 54490E  + 00 

Ms 

3. 

S = 

.90000E+01 

Sl  = 

.49041E+01 

M = 

5. 

S = 

•25000E+02 

Sl  = 

• 13622E+02 

Ms 

10. 

s= 

• 10000E  + 03 

S 1 = 

.54490E+02 

M = 

25. 

s= 

• 62500E+03 

Sl  = 

•34056E+03 

Ms 

100. 

s= 

•10000E+05 

Sl  = 

.54490E+04 

Ms 

250. 

s= 

• 62500E+05 

Sl  = 

• 3405t>E*05 

M = 

1000. 

s= 

. 10000E  + 0 7 

Sl  = 

.54490E+06 

N = sample  size  in  fatigue  test;  n in  the  text 
A (1 ) = shape  parameter  a in  Weibull  distribution 
DZ  = increment  size  dz  for  numerical  integration 
R = fleet  reliability  level 

E (R" ) = expected  value  of  R' ; E[R']  in  Volume  I 
E(R")/R  = reliability  ratio;  E[R']/R  in  Volume  I 
VAR ( R" ) = variance  of  R' 

STD  = standard  deviation  of  R' ; crR, 

VR"  = coefficient  of  variation;  VR|  in  Volume  I 
M = aircraft  fleet  size 

S = scatter  factor,  S,  based  on  R;  from  Eq.  6 
SI  = scatter  factor,  S,  based  on  E [R1 ] ; from  Eq.  6 
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The  output  such  as  that  shown  in  Table  1 was  generated  for 
N = 1,2,...,  10;  A = 0.5,  1,  2,  3,  4,  5,  10;  R = 0.5,  0.6, 

0.7,  0.8,  0.9,  0.99,  0.995,  0.999,  0.9999. 

SECTION  III 

SCATTER  FACTOR  WITH  UNKNOWN  SHAPE  AND  SCALE  PARAMETERS 

(Program  A7720) 

A simplified  flow  chart  of  Program  A7720  is  shown  in 
Figure  2. 

3.1  INPUT  DATA 

There  are  no  input  data  cards  used  for  this  program. 

The  parameters  are  specified  at  the  beginning  of  the  program. 
The  input  parameters  are  as  follows:  R = fleet  reliability, 

M = fleet  size,  NAM  = classification  ASCII  code  for  results, 

N = sample  size  in  fatigue  test;  n in  the  text,  NSMPL  = 
sample  size  in  Monte  Carlo  simulation;  N in  the  text,  10  = 
input/output  peripheral  device  number,  IX  = internal  parameter 
in  RANDU  subroutine  , NST  = a lower  bound  in  array  N,  NND  = 
an  upper  bound  number  in  array  N,  NU  = number  of  divisions  on 
u,  NV0  = number  of  divisions  on  v^,  UU  = upper  cut-off  value 
on  u,  UV0  = upper  cut-off  value  on  Vq . The  parameters  u and 
Vq  are  described  in  Volume  I. 

3.2  PROGRAM  DESCRIPTION 

Program  A7720  consists  of  a main  program  and  six 
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SUBROUTINES 


START 

*4 

INPUT  DATA 


VALUES  FOR  v* , u* , v_ , u 


I 


DISTRIBUTION  OF  V* 


I 


DISTRIBUTION  OF  u* 


DISTRIBUTION  OF  vQ 


DISTRIBUTION  OF  u + 


JOINT  VALUES  OF  u AND  v 
OUTSIDE  SPECIFIED  RANGE 


JOINT  FREQUENCY  F(u,v  ) 

♦ 0 

VALUES  FOR  w,  Z 


0 


DISTRIBUTION  OF  Z 


SAMPL 


DGEN 


ANS 


ANS 


ANS 


ANS 


• ASCEN 
- LIST 
■ ASCEN 
LIST 
ASCEN 
LIST 
ASCEN 


LIST 


l_ 


ANS 


RANDU 

ASCEN 

■* LIST 


E [R"  ] , VAR  R"  , ctr„,  E[R"]/R, 

I 


VR„,  Q* 


END 


Figure  2.  Simplified  Flow  Chart  of  Program  A7720 


1 

I 

I 

I 

RANDU  I 
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subroutines:  SAMPL,  ANS,  RANDU,  ASCEN,  LIST  and  DGEN . The 
program  listing  is  given  in  Appendix  B. 

3.2.1  The  Main  Program  Unit  A7720 
The  main  program  unit  calculates  distributions  of  u, 

Vq,  u*,  vg  and  Z where  these  parameters  are  defined  by  Eqs. 

32,  33  and  36  respectively  in  Volume  I.  The  two 

dimensional  frequencies  of  u and  Vq  are  obtained  by  utilizing 
the  Monte  Carlo  simulation  procedure.  The  reliability 
ratio  E[R"]/R,  the  variance  of  R" , VAR  = E[(R")^]  - (E[R"])^» 
the  standard  deviation  of  R" , aRll  = (VAR)  2,  and  the  coefficient 
of  variation  VR„  = cr  ,,/E  [R"]  are  calculated  from  Eq.  41. 

3.2.2  Subroutine  SAMPL  (U,  VO,  VOS,  IX,  IY,  N,  NSMPL,  IQ) 
Subroutine  SAMPL  is  used  to  generate  values  for  u,  Vq  and 
v^  by  using  the  Maximum  Likelihood  Estimation  (MLE)  procedure. 
The  arguments  for  this  subroutine  are 
U = u,  output  parameter 
VO  = Vq,  output  parameter 
VOS  = v*,  output  parameter 

IX  = any  odd  integer  less  than  nine  digits 
IY  = internal  parameter 

N = sample  size  in  fatigue  test;  n in  Volume  I 
NSMPL  = sample  size  in  Monte  Carlo  simulation; 

N in  Volume  I 

10  = input/output  peripheral  device  number 
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3.2.3  Subroutine  ANS  (A,  NSMPL,  NS,  IQ , NAM1 , NAM2 , M) 


Subroutine  ANS  is  used  to  arrange  the  input  data  in 
proper  order  and  to  print  out  the  probability  distribution. 
The  arguments  for  this  subroutine  are 
A = input/output  array 

NSMPL  = sample  size  in  Monte  Carlo  simulation 
NS  = sample  size 

10  = input/output  peripheral  number 

NAM1  and  NAM2  = classification  parameters  for  the 

results 

M = internal  program  parameter,  or  fleet  size 

3.2.4  Subroutine  RANDU  (IX,  IY,  RAN) 

Subroutine  RANDU  generates  independent  random  numbers 

uniformly  distributed  between  0 and  1.  An  odd  integer, 

IX,  needs  to  be  specified  in  the  main  program.  By  changing 
IX,  a new  sequence  of  random  numbers  is  generated.  The 
arguments  for  this  subroutine  are 

IX  = any  odd  integer  less  than  9 digits 

IY  = internal  parameter  in  the  program 

RAN  = output  uniform  random  number  between  0 and  1. 

3.2.5  Subroutine  ASCEN  (S,  NPT) 

Subroutine  ASCEN  is  used  to  rearrange  given  data  in 

ascending  order.  The  arguments  of  this  subroutine  are 

| 
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S = an  arbitrary  input/output  array 

NPT  = total  number  of  values  in  the  array 

3.2.6  Subroutine  LIST  (A,  NPT,  10,  NAM1 , NAM2 ) 
Subroutine  LIST  is  used  for  listing  the  input  as  a 

probability  distribution.  The  arguments  of  this  subroutine 
are 

A = input  data  arranged  in  ascending  order 
NPT  = total  number  of  values  in  array  A 
10  = output  peripheral  device  number 
NAM1  and  NAM2  = classification  parameters  for 

the  list 

3.2.7  Subroutine  DGEN  (Y,  N,  IX,  RAN) 

Subroutine  DGEN  generates  exponential  random  numbers 

by  utilizing  a uniform  distribution  in  subroutine  RANDU. 

The  arguments  of  this  subroutine  are 

Y = the  exponential  random  number 
N = index  indicating  the  array  of  random  numbers 
IX  = defined  in  RANDU 
RAN  = defined  in  RANDU 


3 . 3 PROGRAM  OUTPUT 

Progam  A7720  prints  out  the  output  as  shown  in 
Tables  2-9. 
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Table  2 Values  for  v^* 


DISTRIBUTION  OF  VO* 

FOR 

SAMPLE 

SIZE  Ns 

3 

SAMPLE 

NO.s 

100 

PROBs 

.0500 

VO* 

.6230 1 49E-0  1 

SAMPLE 

NO.s 

200 

PROB  = 

.1000 

VO* 

z 

. 1618284E+00 

sample 

N0.= 

3 0 0 

PROBs 

.1500 

VO* 

r 

•2759731E+00 

sample 

NO.s 

4 00 

PROBs 

.2000 

VO* 

z 

.3760098E+00 

sample 

N0.= 

soo 

PROBs 

.2500 

VO* 

z 

. 455081 3E +00 

sample 

NO.s 

600 

PROBs 

.3000 

VO* 

z 

.5341889E+00 

sample 

NO.s 

700 

PROBs 

.3500 

VO* 

z 

.6099762E+00 

SAMPLE 

NO.s 

BOO 

PROBs 

.4000 

VO* 

z 

. 6964  087E  + 00 

SAMPLE 

N0.= 

900 

PROBs 

.4500 

VO* 

z 

.7882142E+00 

SAMPLE 

NO.s 

1000 

PROBs 

.5000 

VO* 

z 

. 88035 7 9 E+ 00 

SAMPLE 

NO.s 

1100 

PROBs 

.5500 

VO* 

= 

. 9898615 E +00 

SAMPLE 

NO.s 

1200 

PROBs 

.6000 

VO* 

z 

. 1090346E+01 

SAMPLE 

NO.= 

13  0 0 

PROBs 

.6500 

VO* 

z 

. 1214184E  + 01 

SAMPLE 

NO.s 

1400 

PROBs 

.7000 

VO* 

z 

. 1353251E  + 01 

SAMPLE 

NO.s 

1500 

PROBs 

.7500 

VO* 

z 

. 1553896E  + 01 

SAMPLE 

NO.s 

1600 

PROBs 

.8000 

VO* 

z 

. 1838503E  + 0 1 

SAMPLE 

NO.= 

1700 

PROBs 

.8500 

VO* 

z 

. 2362058E  + 01 

sample 

NO.s 

18  00 

PROBs 

.9000 

VO* 

z 

.3 148309 E +01 

SAMPLE 

NO.s 

1900 

PROB  = 

.9500 

VO* 

z 

.5533883E+01 

SAMPLE 

NO.s 

1 960 

PROBs 

.9800 

VO* 

z 

. 1 4 6 2 6 1 5 E + 0 2 

SAMPLE 

NO.s 

I960 

PROBs 

.9900 

VO* 

Z 

.3651 022E+02 

SAMPLE 

NO.= 

1990 

PROBs 

.9950 

VO* 

Z 

.8491 725E+03 

SAMPLE 

NO.s 

1996 

PR0B  = 

.996  0 

VO* 

- 

. 1077416E  + 08 

SAMPLE 

NO.s 

1 99  B 

PROBs 

.999  0 

VO* 

z 

.7 369510 E + l 1 

sample 

NO.s 

1999 

PROBs 

.9995 

VO* 

z 

. 279407 2E + 14 

VO*  = v* 

0 

SAMPLE  NO.  = a positive  integer  such  that  the  plotting 
position  of  a realization  is  given  by 
(SAMPLE  NO . ) / (NSMPL  +1) 

PROB  = plotting  position 
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Table  3 Values  for  u* 


DISTRIBUTION  OF  U*  FOR  SAMPLE  SIZE  Ns  3 


SAMPLE 

NO.s 

100 

PROBs 

. 0500 

u* 

r 

. 1752809E+00 

sample 

NO  .s 

200 

PROBs 

.1000 

u * 

r 

.2507607E+00 

SAMPLE 

NO  .= 

300 

PR0B  = 

.1500 

u * 

s 

.3056521E+00 

sample 

NO.s 

aoo 

PR0B  = 

.2000 

u* 

s 

.3629199E+00 

SAMPLE 

NO.s 

BOO 

PROBs 

.2500 

u* 

s 

.a2i8aa2E+oo 

sample 

NJ  0 . s 

600 

PROBs 

. 3000 

u* 

s 

,a83a875E+00 

SAMPLE 

NO.S 

700 

PROBs 

.3500 

u* 

s 

. 5392869E  + 0 0 

SAMPLE 

NO.s 

800 

PROBs 

. a 0 0 0 

u* 

s 

.58aia78E+00 

sample 

o 

• 

II 

900 

PROBs 

.9500 

u * 

s 

.6297669E+00 

SAMPLE 

NO.S 

1000 

PROBs 

.5000 

u* 

= 

.6737659E+00 

sample 

NO.s 

1 100 

PROBs 

.5500 

u* 

s 

.7317891E+00 

sample 

NO.s 

12  0 0 

PROBs 

.6000 

tJ* 

s 

.7828305E+00 

sample 

NO  . S 

1300 

PROBs 

.6500 

u * 

- 

.8a05a82E+00 

SAMPLE 

NO  . S 

l aoo 

PROBs 

. 7000 

u* 

s 

.9090222E+00 

SAMPLE 

MJ . s 

150  0 

PROBs 

. 7500 

u* 

z 

.9801 127E+00 

SAMPLE 

II 

• 

o 

16  0 0 

PROBs 

.8000 

IJ  * 

s 

. 1 055699E  + 0 1 

SAMPLE 

NO.S 

1700 

PROBs 

.8500 

u * 

s 

. 1 151 1 12E  + 01 

sample 

NO  . = 

1800 

PROP  = 

.9000 

u * 

s 

. 1286268E+01 

sample 

NO.s 

1900 

PROBs 

.9500 

u* 

s 

. 155aa81Et01 

sample 

NO.s 

196  0 

PROBs 

.9800 

u * 

s 

. 1768360E+01 

sample 

N 0 . = 

1 96  0 

PROBs 

.9900 

u * 

s 

. 1902801E  + 01 

sample 

NCj#  = 

1990 

PRQPs 

.9950 

u* 

s 

.2182300E+01 

sample 

N 0 . = 

1 996 

PROPS 

.9980 

u * 

s 

.2a75597E+01 

SAMPLE 

NU  , - 

1996 

PROBs 

.9990 

Li  * 

s 

.2891  182E  +01 

sample 

1999 

P R U B s 

.9995 

u * 

S 

.2940728E  + 0 1 

U*  = 

U* 

SAMPLE  NO.  = a positive  integer  such  that  the  plotting 
position  of  a realization  is  given  by 
(SAMPLE  NO . ) / (NSMPL  + 1) 

PROB  = plotting  position 
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Table  4 Values  for  Vg 


DISTRIBUTION 

OF  VO 

FQP  ! 

sample  : 

SIZE  Ns 

3 

sample 

NO.s 

100 

PPOB  = 

.0500 

VO 

z 

. 2622875E+00 

SAMPLE 

NO.S 

200 

PPOBs 

.1000 

VO 

z 

.3662895E+00 

SAMPLE. 

NO.s 

300 

PPQB  = 

. 1500 

VO 

s 

. 930770 3E +00 

SAMPLE 

N0.= 

poo 

PPOBs 

.2000 

VO 

z 

. 5051359 E +00 

SAMPLE 

NO.s 

500 

PPOBs 

.2500 

VO 

s 

.5799380F+00 

SAMPLE 

NO.s 

600 

PPOBs 

.3000 

VO 

r 

.6378933E+00 

SAMPLE 

N 0 . = 

700 

PPOB  = 

.3500 

VO 

z 

.7069291E+00 

sample 

NO.S 

600 

PPOBs 

.9000 

VO 

r 

.7713822E+00 

sample 

IMIJ  . S 

900 

PPOBs 

.9500 

VO 

z. 

.8326998E+00 

SAMPLE 

NO.S 

1000 

PR0B  = 

.5000 

VO 

s 

.91991 75E+00 

SAMPLE 

NO.S 

1 100 

PPOBs 

.5500 

VO 

s 

.9893186E+00 

sample 

NO.S 

1200 

PPOBs 

.6000 

VO 

s 

. 1 065538E+01 

sample 

NO.S 

1300 

PPOBs 

.6500 

VO 

r 

. 1 138658E+01 

sample 

N (-1  • = 

1 9 0 0 

PPOBs 

. 7000 

VO 

z 

. 1228765E+01 

sample 

n n . c 

1500 

P P 0 B = 

O 

o 

'T> 

. 

VO 

s 

. 1328859E+01 

SAMPLE 

NO.S 

lt>00 

PPOBs 

.6000 

V 0 

z 

. 1 967891E+01 

SAMPLE 

NO.S 

17  00 

PPOBs 

.6500 

VO 

z 

. 1600763E+01 

sample 

NO.S 

1800 

PPOBs 

.9000 

VO 

z 

•1801288 F +01 

SAMPLE 

NO.s 

190  0 

PPOBs 

.9500 

V 0 

z 

.2160  1 33E  + 01 

SAMPLE 

NO.s 

1 960 

PPOBs 

.9800 

VO 

z 

.2592991E+01 

SAMPLE 

NU.s 

1960 

P P 0 B s 

.9900 

VO 

z 

. 2993689 E + 01 

sample 

NO  . S 

1 990 

PPOBs 

.9950 

VO 

z 

. 3307921 E +0 1 

sample 

NO  . s 

1 99o 

PPOBs 

.9960 

VO 

z. 

. 3659089E  + 0 1 

sample 

N 0 . s 

1 996 

PPOBs 

.9990 

VO 

z 

. 3 8999 36 E +01 

sample 

N 0 • = 

1 999 

PPQbs 

.9995 

\JQ 

z 

.9230395E+01 

SAMPLE  NO.  = a positive  integer  such  that  the  plotting 
position  of  a realization  is  given  by 
(SAMPLE  NO . ) / (NSMPL  + 1) 

PROB  = plotting  position 
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Table  5 Values  for  u 


DISTRIBUTION  OF  U 

FOR 

SAMPLE  SI 

It 

f»  = 

5 

SAMPLE 

NO  .= 

100 

P R O B = 

. 0500 

u 

.6433 0 14 E +00 

SAMPLE 

N 0 . = 

POO 

PR08  = 

.1000 

u 

z 

. 7774431 E + 00 

SAMPLE 

NU  .= 

30  0 

P R U 6 = 

.15  0 0 

u 

z 

.8687P49E+00 

SAMPLE 

NO  . = 

400 

P R 0 B = 

.poo  o 

LJ 

z 

.947P84PE+00 

sample 

NU.= 

500 

PRf)B  = 

. ? 5 0 0 

0 

r 

. 1 OP0291E+01 

SAMPLE 

N0.= 

60  0 

PRC)B  = 

. 3 0 0 0 

LJ 

• 

. 1 106168E+01 

SAMPLE 

NO  . = 

700 

PROB  = 

. 35  0 0 

u 

z 

. 1 189700E+01 

SAMPLE 

NO  .= 

BOO 

P R 0 B = 

.4000 

u 

z 

. 1P77416E+01 

sample 

NO  . = 

900 

PROB  = 

. « 5 0 0 

u 

z 

. 13665P3E+01 

sample 

NO  .= 

10  0 0 

PROB  = 

.5000 

u 

z 

.1484196E+01 

SAMPLE 

N0.= 

1 100 

P R O B = 

.5500 

u 

r 

. 1600597E+01 

sample 

NO  .= 

120  0 

P R 0 B = 

. 6000 

u 

z 

. 171 1895E+01 

SAMPLE 

N0.= 

1300 

PRDB  = 

.6500 

u 

5 

. 185430Pt+01 

sample 

N0.= 

1 POO 

P R O H = 

. 7000 

LJ 

z 

.P068305E+01 

SAMPLE 

NO  .= 

1500 

PRUB  = 

. 750  0 

U 

Z 

. ? 3 7 0 5 4 3 E + 0 1 

SAMPLE 

II 

• 

o 

2 

160  0 

PR0B  = 

.600  0 

u 

Z 

.P759250E+01 

SAMPLE 

NO  .= 

1700 

PROHs 

.6500 

u 

Z 

.3P71693E+01 

SAMPLE 

NO  . s 

1600 

PROB  = 

.9000 

u 

z 

.39878668+01 

SAMPLE 

N0.= 

190  0 

PR()B  = 

.9500 

u 

Z 

. 5 / 0 5 1 P 9 1 + 0 1 

SAMPLE 

N0.= 

1960 

P R 0 B = 

.9800 

u 

z 

. 1 0 01000E  + 0P 

SAMPLE 

NO  .= 

1 98  0 

PR0B  = 

.9900 

u 

z 

. 1581  000E  +0P 

SAMPLE 

NO  .= 

1990 

P R 0 B = 

.9950 

u 

z 

.P361 0G0E+0P 

SAMPLE 

N0.= 

1996 

PRUB  = 

.9980 

u 

z 

.44410 OOE+02 

SAMPLE 

NO.= 

1998 

PROB  = 

.9990 

u 

z 

. 79609988+02 

sample 

N0.= 

1999 

PRQB  = 

.9995 

u 

z 

. 7980998 E +02 

U = u 

SAMPLE  NO.  = a positive  integer  such  that  the  plotting 
position  of  a realization  is  given  by 
(SAMPLE  NO . ) / (NSMPL  + 1) 

PROB  = plotting  position 
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Table  6.  Joint  Values  of  u and  vQ  outside  the  Specified 
Ranges  UU  or  UVO 


Ns 

3 

MS^PL=  1999 

IX=  2810  7 

MJ  = 

50 

NV 0=  ?0  UU= 

25.0000  UV0=  5.0000  DU= 

.5000  D V 0 = 

.2500 

OUTLIERS  f I THER 

IN  U ()P 

.*? 

c 

o 

U = 

.4401E+02 

VOs 

,?682E+00 

J = 

89 

1 = 

2 

U = 

.3981E  tO? 

VOs 

. 6955 E tOO 

J = 

80 

1 = 

3 

U = 

.2691Et02 

V0  = 

. 1355Et01 

J = 

53 

1 = 

6 

IJ  = 

.3221E to? 

VOs 

. Ib53tt01 

J = 

65 

1 = 

7 

u= 

.269  IE tO? 

VOs 

. 5826 E tOO 

J- 

53 

1 = 

3 

Us 

.6361Et0? 

V 0 = 

,6995Et00 

J = 

128 

1 = 

3 

u= 

.7981E+02 

V0  = 

. 1 9 7 9 1 1 0 1 

J = 

160 

1 = 

6 

Us 

.7961Et02 

V 0s 

.67O5Et0Q 

J = 

160 

1 = 

3 

N = sample  size  in  fatigue  test;  n in  Volume  I 

NSMPL  = sample  size  in  Monte  Carlo  simulation; 
N in  Volume  I 

IX  = defined  in  RANDU 

I,  J = array  of  numbers 

NU  = number  of  divisions  on  u 

NVO  = number  of  divisions  on  v^ 

UU  = upper  cut-off  for  u 

UVO  = upper  cut-off  for  v^ 

DU  = increment  on  u,  DU  = UU/NU 

DVO  = increment  on  Vg'  DV0=  UVO/NVO 

U = u 


v0  - vo 


Table  7.  Joint  Frequency  of  u and  vQ,  F ( u,  vQ) 


MU, VO)  f-UK  y A p l t SIZE.  fv  = 


u 

AS 

COLUMN , 

AND 

VO  AS 

1 

t) 

3 

1 

1 

1 

30 

1 09 

69 

62 

59 

96 

22 

69 

9b 

10? 

9 0 

61 

1 1 

'Jb 

66 

59 

57 

33 

5 

25 

2) 

29 

39 

1 H 

5 

1 .1 

20 

Q 

22 

12 

3 

li 

27 

1 9 

12 

9 

3 

6 

) 1 

12 

5 

6 

1 

u 

1 0 

1 0 

6 

3 

2 

1 1 

3 

a 

y 

1 

1 

1 

1 

3 

3 

3 

0 

1 

2 

5 

3 

2 

1 

2 

3 

3 

1 

2 

0 

1 

2 

2 

3 

1 

0 

0 

2 

2 

1 

0 

0 

0 

1 

2 

1 

0 

0 

1 

1 

2 

1 

0 

0 

0 

1 

1 

1 

0 

0 

0 

1 

0 

2 

1 

o 

0 

1 

0 

0 

0 

0 

0 

1 

1 

0 

0 

1 

0 

0 

2 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

1 

1 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

1 

1 

0 

0 

0 

1 

1 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1 

0 

0 

0 

0 

0 

1 

0 

0 

3 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

ROw. 

0 0 0 0 0 0 0 

1*16  5 a 2 0 ? 

9 9 17  17  19  9 ()  ^ 

33  21  10  11  6 3 1 

1 * 1 * 3 6 3 0 1 

U » o 3 1 l i 

13  5 9 l 000 

3 3 13  10  1 

3 2 2 0 0 1 0 

1 0 ) 0 0 0 0 

110  2 0 0 0 

1 o o 0 0 0 

10  10  10  0 

1 0 0 0 0 0 0 

1 0 0 0 0 0 0 

0 1 0 0 0 0 0 

0 0 0 0 0 0 0 

0 0 0 0 0 0 0 

0 0 0 0 0 1 0 

0 o 0 0 0 o 0 

0 1 0 0 0 0 0 

0 0 0 0 0 0 0 

0 0 1 0 0 0 0 

0 0 0 0 0 0 0 

0 0 0 0 0 0 0 

0 0 0 0 0 0 0 

0 o 0 o 0 o 0 

ooooooo 

0 0 0 0 0 0 0 

0 0 0 0 0 0 0 

0 0 0 0 0 0 0 

ooooooo 

0 0 0 0 0 0 0 

OOOOOOO 
0 0 0 0 0 0 0 

o o 0 o 0 0 0 

OOOOOOO 
0 0 0 0 0 0 0 

OOOOOOO 
OOOOOOO 


0 

1 

2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


0 0 

0 o 

1 0 

1 0 

0 0 

1 1 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 o 

0 0 

0 0 

0 0 

0 0 

0 o 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

1 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

0 0 

1 0 
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Table  8 Values  for  Z 


DISTRIBUTION-  OF-  Z FOR  FLEET  SJZE  Ms  1 


SAMPLE 

NO  . = 

100 

P k 0 B = 

.0600 

Z 

Z 

-.  1627202E+01 

SAMPLE 

N 0 . - 

200 

PROtJ  = 

.10  0 0 

l 

Z 

-.86816336+00 

SAMPLE 

NO.  a 

3 0 0 

P P t.i  B = 

.16  0 0 

l 

z 

-.60386236+00 

SAMPLE 

N 0 . = 

400 

P P U b = 

.2000 

L 

z 

-.96366068+00 

sample 

NO.  a 

50  0 

P P 0 H a 

.2600 

7 

z 

-.39069698+00 

SAMPLE 

0 u . a 

600 

PPOBa 

. 3000 

L 

z 

-. 2223126 E+00 

Sample 

NO.  a 

7 00 

PPOBa 

.3600 

1 

z 

-.  1006370E+00 

Sample 

NO.  a 

6 0 0 

p p r i r s 

.900  0 

l 

s 

-.  1 083620  6-  01 

SAMPLE 

NO  . = 

900 

P k 0 B a 

.9600 

1 

z 

.92152226-01 

SAMPLE 

NO.  a 

1000 

P P 0 B a 

.6000 

2 

z 

. 1 713892E  + 00 

sample 

NO.  a 

1100 

P P 0 B a 

.6600 

1 

= 

.27988486+00 

SAMPLE 

NO.  a 

1200 

P B 0 h a 

.600  0 

1 

z 

. 3893890E  + 00 

SAMPLE 

Ml.s 

1 300 

PPOBa 

.6600 

2 

z 

.52399206+00 

sample 

NO.  a 

1 A 0 0 

P k Q 6 = 

.7000 

1 

z 

.6791 7286+00 

SAMPLE 

N 0 . S 

1 6 0 0 

P K O H a 

.7600 

2 

z 

.67807206+00 

sample 

NO.  8 

160  0 

P P 0 B a 

. 6 0 0 0 

7 

z 

. 1 l97860t  + 0 1 

sample 

NO.  a 

1700 

PPOBa 

.6600 

7 

z 

.19173876+01 

Sample 

NO.  a 

160  0 

PPOBa 

.9000 

Z 

z 

.18998016+01 

SAMPLE- 

NC'.S 

1900 

PPOBa 

.9600 

2 

z 

.31107526+01 

SAMPLE 

NU.s 

I960 

P k 0 H a 

.960  0 

l 

z 

.96793016+01 

SAMPLE 

NO.  a 

1 9 6 0 

PP0B  = 

.9900 

1 

z 

.61066146+01 

SAMPLE 

i'll-  .= 

1 990 

PPOBa 

.9960 

l 

z 

.77697896+01 

sample 

NO.  a 

1996 

P P 0 P = 

.9960 

7 

z 

.2141 146E+02 

sample 

ivU  . = 

1996 

PPOb  = 

.9990 

L 

z 

.27566676+02 

sample 

NU.s 

1 999 

P P 0 B = 

.9996 

L 

z 

.29619786+02 

Z = Z 

M = aircraft  fleet  size 

SAMPLE  NO.  = a positive  integer  such  that  the  plotting 
position  of  a realization  is  given  by 
(SAMPLE  NO. ) / (NSMPL  +1) 

PROB  = plotting  position 
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Table  9 

Expected 

Value , 

Variance , 

Standard 

Deviation, 

Reliability  Ratio,  and  Coefficient 

of  Variation  of  R" . 

3 !'■ 

8 V H l = J 9 9 9 * = 

1 

K 

E(R""  ) 

UK 

8 T D 

E (R"")/R 

V K "" 

QS 

. b 0 0 (i 

.46  7 6 

.210/ 

. 9 7 b 7 

.9316 

.139476+01 

.600  0 

.b  707 

.0031 

. 2076 

. 9 b 1 2 

.3637 

.21 1366+01 

. 7 0 0 0 

. 6520 

.0362 

.19  0 3 

.931b 

.2918 

.34162E+01 

. 6 0 (i  0 

.7323 

. 0296 

.17  2 0 

.9153 

.2399 

.64077E+01 

.9000 

. 6 J 7 3 

.01  74 

. 1 320 

.9062 

.1616 

. 16770E+02 

. 9 b 0 0 

.6661 

.0106 

.10  9 0 

.9136 

.1196 

. 36765  F.  + 02 

.9900 

. 9 ? 6 0 

.0097 

. 0669 

.9379 

. 0793 

.25341F+03 

. 999  (i 

.963/ 

. 0020 

. 0 9 9 9 

.9697 

. 0465 

.4961 3E +04 

.9999 

. 9 7 8 41 

.0013 

. 0 3 b 6 

.9785 

. 0364 

•34388 E +05 

N = sample  size  in  fatigue  test;  n in  Volume  I 

NSMPL  = sample  size  in  Monte  Carlo  simulation; 

N in  Volume  I 

M = aircraft  fleet  size 

R = assigned  reliability  values 

E(R"")  = expected  values  of  R"  obtained  from  Eq.  41. 

2 2 

VAR  = E[  i(R" ) ] - E [R"],  variance  of  R"  in  Volume  I 

STD  = aRll , standard  deviation  of  R"  in  Volume  I 
E(R"")/R=  reliability  ratio 

VRIII1  = coefficient  of  variation  of  R" , STD/Ef  R"  ] volume  I 
QS  = Q*  in  Volume  I 
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Appendix  A 


PROGRAM  A7701(  OUTPUT  * T APE6  = 0UTPUT ) 

PURPOSE 

SCATTER  FACTOR  CALCULATION 

WITH  KNOWN  WEIBULL  SHAPE  PARAMETER. 

DESCRIPTION  OF  THE  PROGRAM 

THIS  PROGRAM  CALCULATES  E(R»).  E(R»)/R. 

VAR(R»).  VR"  FOR  GIVEN  WEIBULL  SHAPE  PARAMETER  A 
SAMPLE  SIZE  N»  FLEET  SIZE  AM.  AND  RELIABILITY  LEVEL  R. 
IT  ALSO  CALCULATES  CORRESPONDING  SCATTER  FACTORS 
S BASED  ON  R AND  SI  BASED  ON  E(R"). 

REMARKS 

INTEGRATION  TO  GET  E(R">  AND  E((R'')**2>  IS  CARRIED 
OUT  BY  SIMPSON"S  METHOD. 

DATE  7/16/1977 

DIMENSION  X (1000) .A (50) »CV (50) .STD (50) * Z < 1 000 ) ♦ NN ( 50 > 
DIMENSION  XX (1000) *R(50) .ERR (50) * ARG3 ( 1 000 ) . ARG2 ( 1000 ) 
DIMENSION  AM (50) .ANZA(IOOO) .VAR (50) 

ASSIGNMENT  OF  SPECIFIED  ARRAY  INPUT  DATA. 

AM=  FLEET  SIZE 
A=  WEIBULL  SHAPE  PARAMETER 
R=  SPECIFIED  RELIABILITY  LEVEL 
NN=  SAMPLE  SIZE 

DATA  AM/1 . .3. *5. .10. .25. *100. .250. , 1000./ 

DATA  A/0. 5.1. .2. .3. .4. .5.  .10./ 

DATA  R/0. 5. 0.6. 0.7 .0.8. 0.9. 0.99 *0.995 .0,999. 0.9999/ 

DATA  NN/ 1 . 2 ♦ 3 . A . 5 . 6 * 7 . 8 * 9 . 1 0 . 20/ 

ASSIGNMENT  OF  ARRAY  SIZE  TO  BE  USED  (UP  TO  50). 

IN  — FOR  SAMPLE  SIZE  NN 
I M — FOR  FLEET  SIZE  AM 
IA  — FOR  WEIBULL  SHAPE  A 
IR  — FOR  FLEET  RELIABILITY 

IN=  1 1 

I M = 8 
IA  = 7 
I R = 9 


TO  SPECIFY  OUTPUT  PERIPHERAL  NUMBER  10. 

10  = 6 

TO  SPECIFY  NUMBER  OF  INCREMENTS  NPT  FOR  SIMPS0N"S 
NUMERICAL  INTEGRATION  (UP  TO  1000). 

USE  NPT=400  IN  THIS  PROGRAM. 

NPT  =400 
NPT1=NPT-1 

TO  CALCULATE  FOR  EACH  SAMPLE  SIZE  N. 

2B0  DO  200  JN= 1 . 1 N 
N=NN ( JN) 

AN=FLOAT (N) 

ANl=l ./AN 

N2=2*N 

ANN=AN**N 

CALL  GAMMA  (N.6A) 
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TO  calculate  FOR  EACH  w e i bull  shape  parameter  a. 

DO  90  K = 1 » I A 
DZ  = 0 • 025 

IF  (K  .GE.  3)  DZ  = 0 ♦ 0 1 
PRINT  25,  N , K * A ( K ) ,DZ 

25  FORMAT ( 1H1/////5X,*  SAMPLE  SIZE  N=*,I5// 

15X,*  WE  I BULL  SHAPE  A (», 13,*) =*»F10.5» 

25X,»  INCREMENT  ON  Z DZ=*,F10.5/) 

A 1 = 1 ./A (K) 

AKN1=A (K) »ANN/GA 
ATN1=A (K) *AN-1 . 

TO  CALCULATE  FOR  EACH  RELIABILITY  LEVEL  R. 

DO  50  L= 1 * I R 
RR= ALOG ( R ( L ) ) 

PRINT  110,  R ( L ) 

110  FORMAT (/5X»*  FLEET  RELIABILITY  LEVEL  R=*,F10.5) 

TO  CALCULATE  ARGUMENTS  IN  EOS. (22)  AND  (23) 

BY  USING  DENSITY  FUNCTION  OF  Z. 

DO  60  1 = 1, NPT 
IF  < L ,GT.  1)  GO  TO  AO 
Z(I)=FLOAT(I)*DZ 
ARG2 ( I ) =-AN*Z ( I ) **A ( K ) 

IF  ( ARG2 ( I ) ,LT.  -100.)  GO  TO  AO 
ARG3 ( I ) =AKN1*Z ( I ) tt*ATNl°EXP(ARG2(I ) ) 
ANZA(I)=AN*Z(I)**A(K) 

AO  IF  ( ARG2 ( I ) .LT.  -100.)  GO  TO  30 
ARG=1./(1.-RR/ANZA(I) ) 

80  X ( I ) =ARG*»N*ARG3 ( I ) 

XX ( I ) =X ( I ) *ARG**N 
GO  TO  60 
30  X { I ) =0 . 

XX ( I ) =0. 

60  CONTINUE 

C TO  PERFORM  INTEGRATION  TO  ESTIMATE  E ( R" ) AND  E((R»)**2). 

CALL  SIMPN(X*NPT1,0Z,AREA) 

CALL  SIMPN(XX»NPT1»DZ»AREX) 

C 

C TO  CALCULATE  E(R")/R,  VAR(R"),  STD ( =STANDARD  DEVIATION), 

C AND  VR" (=COEFFICIENT  OF  VARIATION)  AND  ALSO  PRINT  OUT. 

ERR ( L ) = ARE A/R ( L ) 

VAR(L) =AREX-AREA»*2 
IF (VAR(L) .LT.O.)  VAR(L)=0. 

STD(L)=SQRT (VAR(L) ) 

CV(L)=STD(L)/AREA 

PRINT  10,  AREA,ERR(L) , V AR ( L ) »STD(L) ,CV(L) 

10  FORMAT ( 9X * *E ( R" ) = * ♦ E 1 2. 5 ♦ 5X , *E ( R" ) /R  = * , E 1 2. 5/ 

19X»»VAR(R")=»,E12.5,5X,»STD  =*,E12.5/ 

29X»»VR"  =»»E12.5/ 

318X»*FLEET  SI ZE» ,5X ♦ *S. F . BASED  ON  R*,AX, 

A*S • F • BASED  ON  E(R'»)*) 

C 

C TO  CALCULATE  SCATTER  FACTORS  S BASED  ON  R 

C AND  SI  BASED  ON  E(R»>  FOR  EACH  FLEET  SIZE. 
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DO  70  J= 1 * I M 

S1=(AM(J)/(AN*(1./AREA*»AN1-1.) ) )**A1 
SS=<AM(J)/<ANM1./R<L>*»AN1-1.) ) >**A1 
PRINT  120  » AM ( J) *SS»S1 

120  FORMAT ( 17X»*  M=*,F6.0*7X**S=  **E12.5*4X**S1=  **E12.5) 
70  CONTINUE 
50  CONTINUE 
90  CONTINUE 
200  CONTINUE 
STOP 
END 

SUBROUTINE  GAMMA (N*GA) 

TO  CALCULATE  GAMMA  FUNCTION. 

USE  FACTORIAL  RELATIONSHIP  SINCE  N IS 
AN  INTEGER  IN  EVERY  CASE. 

N MUST  BE  LESS  THAN  50. 

GA= 1.0 

IF(N.EO.l)  RETURN 
N1=N-1 

DO  10  1 = 1 * Nl 
10  6A=GA*FL0AT ( I ) 

RETURN 

END 


SUBROUTINE  SIMPN  < AF , NPOI NT  * DSTEP ♦ ARE A ) 
DIMENSION  AF(1) 

TO  CARRY  OUT  INTEGRATION  FOR  GIVEN  VALUES  OF 
FUNCTION  AF. 

NP=(NP0INT-l)/2 

MN=NP-1 

ODD  = 0. 

EVEN=0. 

END=AF(l)*AF(NPOINT) 

DO  10  I=1*NP 
10  ODD=ODD+ AF (2*1) 

DO  20  1 = 1* MN 
20  EVEN=EVEN*AF (2*1*1) 

AREA=(4.0*ODD+2.0*EVEN*END)  »DSTEP/3.0 

RETURN 

END 
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Appendix  B 


PROGRAM  A7720  (OUTPUT*  T APE 6=0UTPUT ) 

C PUPPOSE 

C SCATTER  FACTOR  CALCULATION 

C WITH  UNKNOWN  WE  I BULL  SHAPE  AND  SCALE  PARAMETERS. 

C DESCRIPTION  OF  THE  PROGRAM 

c this  program  first  generates  nsmpl  number  of 

C SAMPLE  sets  of  size  n for  u,  U*,  VO  AND  vo*» 

C AND  THEN  CALCULATES  THEIR  PROBABILITY  DISTRIBUTION 

C AS  WELL  AS  JOINT  DENSITY  F(U»V0). 

C SECOND,  IT  CALCULATES  PROBABILITY  DISTRIBUTION 

C OF  Z FOR  SAMPLE  SIZE  N AND  FLEET  SIZE  M BY  GENERATING 

C CORRESPONDING  SAMPLES  OF  W. 

C IT  FINALLY  CALCULATES  E ( R " ) , VAR (RM ) , VR".  ECR'M/R  AND  U*. 

C remarks 

C IN  THE  PROGRAM,  »*»  IS  REPLACED  BY  "S"»  EXAMPLE,  U#=US. 

C DATE  10/25/77 

C 

DIMENSION  0(1999) ,V0<1999) , VOS (1999) ,Z(1999) 

DIMENSION  NVOU (20,50), K(9),M(8),N(10>  ,NAM(6) 

EQUIVALENCE  (VOS,Z) 

C 

C ASSIGNMENT  OF  SPECIFIED  ARRAY  INPUT  DATA. 

C N=  SAMPLE  SIZE. 

C M=  FLEET  SIZE. 

C R=  SPECIFIED  RELIABILITY  lEVEL. 

C NAM=  DATA  CLASSIFICATION  ASCII  CODE. 

DATA  N/2,3,4,5,6, 7,8,9,10,20/ 

DATA  M/ 1,3,5, 10,25, 100,250, 1000/ 

DATA  R/0.5,0. 6,0.7, 0.8,0.9,0.95, 0.99, 0.999, 0.9999/ 

DATA  NAM/2HU  ,2H  * 2HU* » 2HV 0 , 2H*  ,2HZ  / 

C 

C INPUT  DATA  ASSIGNMENT. 

C NSMPL=  NUMBER  OF  SAMPLE  SETS  OF  SIZE  N. 

C I X=  INITIAL  VALUES  FOR  UNIFORM  RANDOM  NUMBER 

C GENERATION  SUBROUTINE  RANDU ( I X , I Y , RAN) . 

C NST  AND  NND=  CONTROL  PARAMETER  FOR  DO-LOOP  CALCULATION 

C FOR  EACH  SAMPLE  SIZE  N. 

C NST=  LOWER  BOUND  NUMBER  IN  ARRAY  N. 

C NND=  UPPER  BOUND  NUMBER  IN  ARRAY  N. 

C NU=  ARRAY  SIZE  FOR  U. 

C NV0=  ARRAY  SIZE  FOR  VO. 

C UU=  UPPER  CUT-OFF  VALJE  FOR  U. 

C U V 0 = UPPER  CUT-OFF  VA.UE  FOR  VO. 

C LU  AND  10=  SELECT  CODES  OF  DATA  INPUT/OUTPUT  PERIPHERAL. 

C 

NSMPL= 1 999 
IX=23451 
NST  = 1 

nnd=io 

NU  = 50 
NV0=20 
UU=25 • 

UV0=5. 

CALCULATION  FOR  EACH  SAMPLE  SIZE  N ( KN ) 

FNSMP1=FL0AT ( NSMPL* 1) 

FNSMP?=1 .O/FNSMPl 
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DO  100  KN=NST  , NND 
NS=M(KN) 

C 

C TO  GENERATE  NSMPL  NUMBER  OF  SAMPLE  SETS  OF  SIZE  N 

C FOR  U,  VO,  AND  VO*. 

CALL  SAMPL (U. VO.VOSt IX.NS.NSMOLt 10) 

C 

C TO  CALCULATE  PROBABILITY  DISTRIBUTION  OF  VO* 

C BY  ARRANGING  IN  ASCENDING  ORDER. 

CALL  ANS(V0S*NSMPL*NS»I0»NAM(4) ,NAM<5) ,0) 

C 

C TO  CALCULATE  PROBABILITY  DISTRIBUTION  OF  U*. 

C BECAUSE  OF  LlMlTTED  MEMORY  SIZE  OF  THE  COMPUTER  USED* 

C AND  BECAUSE  OF  VALUES  U AND  VO  TO  BE  RESERVED  FOR 

C CALCULATION  OF  Z.  ARRAY  VOS  IS  USED  AS  COMMON  ARRAY. 

DO  110  1=1, NSMPL 
110  V0S(I)=1./U(I) 

CALL  ANSI  VOS, NSMPL, NS, 10, NAM (3) »NAM(2) ,0) 

C 

C TO  calculate  probability  distribution  OF  VO. 

DO  120  1=1, NSMPL 
120  V0S(I)=V0(I) 

CALL  ANS(VOS»NSMPL*NS»IO»NAM(4) ,NAM(2) ,0) 

C 

C TO  CALCULATE  PROBABILITY  DISTRIBUTION  OF  U. 

DO  130  1=1, NSMPL 
130  V0S(I)=U(I) 

CALL  ANS (VOS  , NSMPL , NS , I 0 , N AM ( 1 ) , NAM ( 2 ) ♦ 0 ) 

C 

C TO  CALCULATE  JOINT  DENSITY  F(U,V0),  E ( R" ) , VA«(K"), 

C VR",  E(R")/R»  AND  Q*. 

C 

DU=  UU/FLO A 1 (NU) 

DVO=UVO/FLOAT (NVO) 

DO  150  J= 1 ♦ NU 
DO  150  1 = 1, NVO 
NVOU ( I . J) =0 
150  CONTINUE 
C 

C TO  PRINT  OUT  NECESSARY  INFORMATIONS  FOR  CHECK  PURPOSE 

WRI TE (6,600)  NS, NSMPL, IX 

600  FORMAT (1H1/////5X,*  N=*,I5,*  NSMPL=*»I5,*  IX=*»I20) 
WR I TE ( 6 » 6 1 0 ) NU,NV0,Ul),JV0,DU,DV0 

610  format  t /sx , * nu=*,I5,*  nvo  = *,I5,*  uu=*,fb.a.*  uvo=*, 
lFB.A,*  DU=*»F8.A»*  DV0  = *,FB.  A // ) 

WRI TE (6,61 1 ) 

611  FORMAT ( /5X » * OUTLIERS  EITHER  IN  U OR  IN  VO*/) 

DO  160  K=l» NSMPL 

J=IFIX (U(K) /DU) +1 
I = IFIX ( VO  <K) /DVO) ♦ 1 

IF  ( J.LE.NU  .AND.  I.lE.NVO)  GO  TO  10 
C 

C TO  PRINT  OUT  OUTLIERS  EITHER  IN  U OR  VO 

C FROM  THE  PRE-DEFINED  RANGE. 

vRITE (6,630)  U(K) , VO (*) , J, I 

630  FORMAT*  5X , * U=* , E 1 2 . A ♦ 3X , * VO =* , E 1 2 . A , 3X , * J=*,I6,3X, 
1*  I=*. Ih) 
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IF(J.GT.NU)  J = NJ 
IF(I.GT.NVO)  1 =NV0 
10  NVOU ( I » J) =NV0U ( 1 » J) +1 
160  CONTINUE 

TO  PRINT  OUT  JOINT  FREQUENCY  F(U,V0). 

JOINT  FREQUENCY  F(U»V0)  IS  LISTED  WITH  IJ  AS  COLUMN, 
AND  VO  AS  ROW. 

WRITE(6,640)  NS 

640  FORMAT ( 1H1/////5X,*  F(U,V0)  FOR  SAMPLE  SIZE  N=*»I5/ 

1 8X » * U AS  COLUMN,  AND  VO  AS  ROW.*/) 

WRITE ( 6 , 650 ) NVOU 
650  FORMAT  <20  < IX, 13) ) 

FOR  BRIEF  VERSION  OF  LISTING  OF  F(U,V0). 

WRITE  < 6 , 64  0 ) NS 
DO  170  J = 1,40 

WR I TE ( 6 , 660 ) < NVOU ( I , J ) , I = 1 , 1 6 ) 

660  FORMAT (6X»16(1X,I3) ) 

170  CONTINUE 

TO  CALCULATE  DISTRIBUTION  OF  Z FOR  EACH  FLEET  SIZE. 
SAME  INITIAL  VALUES  FOR  UNIFORM  RANDOM  NUMBER 
GENERATION  SHALL  BE  USE 0 FOR  EACH  CASE  OF  FLEET  SIZE. 
IXS=IX 

DO  200  KK=1,8 
MM=M ( KK ) 

E!M  = FLOAT  (MM) 

ENS=FLOAT (.NSMPL) *1. 

IX=IXS 

DO  210  1=1, NSMPL 
CALL  RANOU  (IX, IT, HAN) 

I X = I Y 

W=-ALOG (RAN) /FM 
V0W=V0(I)/W 

IF  (VOW  . LT . l.E-10)  VO  W = 1 . E- 1 0 
Z ( I ) =U ( I ) *ALOG10 (VOW) 

210  CONTINUE 

CALCULATE  DISTRIBUTION  OF  Z FOR  EACH  FLEET  SIZE 
CALL  ANS(Z,NSMPL,NS,I0,NAM(6) ,NAM(2) , MM ) 

TO  DEFINE  THE  CASE  CALCULATED. 

WRITE  (6,670)  NS, NSMPL, MM 

670  FORMAT < 1H1///5X,*  N=*,I5,*  NSMPL=*»I6,*  M=*,I6> 
WRITE (6,680) 

680  FORMAT ( //8X » * R*,7X,*  E ( R"* ,*")*, 3X , * VAR*,6X, 

1*  STD* , 3X » * E (R"*,*") /R»*  IX,*  VR"* ,*••*, 7X ♦ * QS*/) 

TO  ESTIMATE  E(R"),  VAR ( R") , VR",  E(R")/R,  AND  0*. 

DO  300  K=1,R 
FMALR=FM/ALOG ( 1 ,/R (K) ) 

ER  = 0 . 0 
VAR  = 0 » 0 
DU2=DU/2 . 

DV02=DV0/2. 

DO  310  J=1,NU 
DO  320  1=1, NVO 


25 


o n 


IF  (NVOU  ( I « J)  .ta.  0)  GD  TO  320 
lJX  = FLOAT  < J)*DU-0U2 

V 0 X = F L 0 A r ( I ) *UV0-DV02 
PARM=VOX*FMALR 
20=UX*AL0610 (FARM) 

PAR2=FL0AT  ( iM V TMJ  ( I ♦ J ) >/FMS 
DO  330  L = 1 « NSMPL 

IF  (ZO  .L£ . Z (L) ) GO  TO  20 
330  CONTINUE 
20  FZ  = FLOAT  (U/FMS 
F Z2  = F Z**2 
ER=ER*FZ*PAR2 

V AR  = VARt-FZ2*PAR2 

320  continue 
310  CONTINUE 

VAR=VAR-ER"*2 
STD=SQRT (VAR) 

Cy/  = STD/ER 
RATI 0=£R/R ( K ) 

ZER=ER*F  NSMP1 
IZ=IFIX (ZER) 

IZ1=IZ*I 

FIZ=FLOAT ( 1Z) /FNSMP1 

ZANS=Z ( IZ) ♦ (ER-FIZ) * ( Z ( IZ1 ) -Z ( 12) ) /FN5MP2 
QS=10.»*ZANS 

WR I TE ( 6 ♦ by 0 ) R(<)*£R»VAR*STD* RATIO* CV»QS 
6R0  FORMAT  (6X»&F9.<+,  1X*E12.  o) 

300  CONTINUE 
200  CONTINUE 
100  CONTINUE 
STOP 
END 

SUBROUTINE  ANS ( A » NSMPL * NS ♦ I D * N AM  1 « NAM2 ♦ M ) 
DIMENSION  A ( 1 ) 

TO  aR I NT  OUT  PROBABILITT  DISTRIBUTION 
FOR  GIVEN  ARRAY  A. 

CALL  ASCEN  < A * NSMPL ) 

IF  (M  .EO.  0)  GO  TO  10 
WRITE (6*100)  NAM1 *NAM2*  M 
100  FORMAT ( 1H1/////SX**  DISTRIBUTION  OF  *»2A2* 
1*  FDR  FLEET  SIZE  M=»*I6//) 

GO  TO  20 

10  WRI TE ( 6 * 1 1 0 ) NAM1 »NAM2»  NS 
110  FORMAT (1H1/////5X**  DISTRIBUTION  OF  **2A2* 
1*  F DR  SAMPLE  SIZE  N=**I3//) 

20  CALL  LIST  (A.NSMPL* I0*NAM1 , NAM2) 

RETURN 

END 
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SUBROUTINE  L I ST  ( A ♦ NPT  ♦ I 0 ♦ NAM1  * NAM2 ) 

DIMENSION  A ( 1 ) *P (6) 

TO  LIST  CORRESPONDING  REALIZATION  VALUES  FOR 
SPECIFIC  PROBABILITY  LEVELS. 

DATA  P/0. 98. 0. 99 * 0. *96* 0.998*0. 999 *0.9995/ 

NPT1 =NPT ♦ 1 
ANPTl=FLOAT (NPT1) 

ND= ( NPT ♦ 1 ) /20 
IF  (ND  .LT.  1)  ND= 1 
DO  10  I =ND ♦ NPT . NO 
PROB=FLOAT (IJ/ANPT1 
WRITE (6. 100)  I*PROB»NAMl»NAM2*A(I) 

10  CONTINUE 
DO  20  1=1.6 
J=IFIX(ANPT1*P(I)  ) 

WRirE(6*100)  J.P(I)  . NAM1 » NA  M2  » A ( J ) 

IFIJ.EU.NPT)  GO  TO  30 
20  CONTINUE 

100  FORMAT ( 5X ♦ * SAMPLE  NO .=* . I 5 » 3X ♦ * PROB=**F7.4.3X» 
12A2 . * =* ♦ E 1 4 • 7 ) 

30  RETURN 
END 

SUBROUTINE  ASCEN(S.NPT) 

DIMENSION  S ( 1 ) 

TO  ARRANGE  GIVEN  ARRAY  5 IN  ASCENDING  ORDER 
IF  (NPT  .LT.  2)  RETURN 
N2=NPT/2 
DO  20  J= 1 ♦ N2 
m*i=NPT-J+1 
K = J 
M = MM 

DO  10  I = J ♦ MM 
IF  (S(I)  .LT.  S ( K ) ) K = I 
IF  (S(I)  .GT.  S ( M ) ) M=I 
10  CONTINUE 
C TEST  FOR  J=K 

IF  (J  .NE.  K)  GO  TO  50 
IF  (MM  .EU.  M)  GO  TO  20 
TEMP=S(M) 

S (M) =S(MM) 

S (MM) =TEMP 
GO  TO  20 

C TEST  FOR  MM=M 

50  IF  (MM  .NE.  M)  GO  TO  60 

TEMP=S(J) 

S ( J) =S (K) 

S(K) =TEMP 
GO  TO  20 

60  IF  (J  .EQ.  M)  60  TO  30 

C J/=M.MM/=K 

TEMP=S( J) 
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S(J)=S(K) 

S ( K ) =TEMP 
TEM?1=S (MM) 

S(MM)=S(M) 

S ( M ) = TEMP  1 
GO  TO  20 

30  IE (MM  . EQ.  K)  GO  TO  40 
C J = M 

TEMa 1 =S ( MM ) 

S(MM)=S(M) 

S ( M ) =TEMP 1 
TEMP=S(J) 

S(J)=S(K) 

S (K) =TEMP 
GO  TO  20 
C J=M»MM=K 

40  TEMP=S(J) 

S(J)=S(MM) 

S (MM) =TEMP 
dO  CONTINUE 
PE  TURN 
END 


SUBROUTINE  S AMPL ( U , V 0 * V 0 S * I X , N , NSMPL » I O ) 

DIMENSION  IJ  ( 1 ) , V0  ( 1 ) , VOS<  1 ) ♦ Y (20)  ♦ RAN  (20)  * ALY  (20) 

TO  GENERATE  NSMPL  NUMBER  OF  SAMPLE  SETS  OF  SIZE  N 
FOR  U»  VO  AND  VO*. 

SAMPLE  SIZE  N MUST  BE  LESS  THAN  OR  EQUAL  TO  20  FOR 
MINI-COMPUTER. 

DU  = 0 . 2 
NPT  = 400 

DO  100  L=  1 » NSMPL 
CALL  DGEN(Y»N»IX»RAN) 

KASE=0 

PAR5=0. 

00  110  1 = 1, N 
ALY ( I ) = ALOG  ( Y ( I ) ) 

PAR5=PAR5+ ALY  ( I ) 

110  CONTINUE 

PAR5=PAR5/FL0AT (N) 

UM I N— 0 • 0 1 

FIND  THE  SOLUTION  OF  PHl(U)  IN  MLE  EQUATION. 

00  120  J= 1 , NPT 
Ul=FLOAT  < J-l > *DJ*UMIN 
10  PAR  1 = 0. 

PAR2=0 . 

DO  130  1 = 1,  N 
Y I U = Y ( I ) **U 1 
PAR1=PAR1+YIU*ALY ( I ) 

PAR2=PAR2*YIU 
130  CONTINUE 

PHI2=PAR1/PAR2-1 ./U1-PAR5 

IF  ( ABS (PHI?) .LT.  0.001)  GO  TO  20 

IF  (J  .EQ.  1)  GO  TO  30 

IF  (KASE  .EQ.  1)  GO  TO  20 

IF  ( PH 1 1 .GT.  0.)  GO  TO  40 

IF  (PHI 2 .LT.  0.)  GO  TO  30 
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GO  TO  50 

*f0  IF  (PhI2  .GT.  0.)  GO  TO  30 
50  Ul  = J1 -PHI  2 / ( PHI  2-PHI  1 ) # i)U 
KASE=1 
GO  TO  10 

30  PHI1=PHI2 
120  CONTINUE 
20  U ( L ) = U 1 
US=1./U1 

VOS(L)=PAH2/FLOAT(N) 

UV0=US*ALOG10 (V0S(L) ) 

IF  (UVO  .GT.  -30.)  GO  TO  60 
VO (L) =0. 

GO  TO  100 

60  IF  (UVO  .LT.  30.)  GO  TO  70 
VO (L>  =1 .E30 
GO  TO  100 

70  VO (L) =V0S (L) **US 

100  CONTINUE 

RETURN 

END 

SUBROUTINE  DGEN ( Y • N ♦ I X » R AN ) 

DIMENSION  Y ( 1 ) ♦ RAN ( 1 ) 

SUBROUTINE  TO  GENERATE  EXPONENTIAL  RANDOM  NUMBER 
OF  SIZE  N 

BY  Y = -LN  ( X ) ♦ inlHERE  X --  UNIFORM  RANDOM  NUMBER 

DO  10  1=1 »N 

CALL  RANDU ( IX  * I Y »RAN ( I ) ) 

I X = I Y 

10  Y (I)=-ALOG(RAN(I) ) 

RETURN 

END 


SUBROUTINE  RANDU ( IX * I Y » RAN) 

M=2**32 

K=2**7*l 

IY=MOD(IX*K*l«M> 

RAN=FLOAT ( 1 Y ) /FLOAT ( M ) 

RETURN 

END 
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